Van der Waals Density Functional for General Geometries 
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A scheme within density functional theory is proposed that provides a practical way to generalize 
to unrestricted geometries the method applied with some success to layered geometries [H. Rydberg, 
et at, Phys. Rev. Lett. 91, 126402 (2003)]. It includes van der Waals forces in a seamless fashion. 
By expansion to second order in a carefully chosen quantity contained in the long range part of 
the correlation functional, the nonlocal correlations are expressed in terms of a density-density 
interaction formula. It contains a relatively simple parametrized kernel, with parameters determined 
by the local density and its gradient. The proposed functional is applied to rare gas and benzene 
dimers, where it is shown to give a realistic description. 

PACS numbers: 31.15.Ew,71.15.Mb,61.50.Lt 



Density functional theory (DFT) for molecules and 
materials is widely applied with approximate local and 
semilocal density functionals for the interaction effects. 
For largely homogeneous systems, for example, simple 
metals and semiconductors, the local-density approxima- 
tion (LDA) for these effects is appropriate. For inho- 
mogeneous systems, for example, transition metals, ionic 
crystals, compound metals, surfaces, interfaces, and some 
chemical systems, semilocal-density approximations, like 
members of the family of generalized gradient approxi- 
mations (GGA), work well. Today DFT describes cohe- 
sion, bonds, structures, and other properties very well for 
dense molecules and materials, as shown by recent stud- 
ies for both single molecules [l| and dense solid-state |2j 
systems. However, sparse systems, including soft mat- 
ter, van der Waals complexes, and biomolecules, are at 
least as abundant. They have interparticle separations, 
for which nonlocal, long-ranged interactions, such as van 
der Waals (vdW) forces, are influential. 

The aim of this letter is to develop and apply a van der 
Waals density functional (vdW-DF) for general geome- 
tries to supplement the planar vdW-DF that we recently 
applied with some success sf to several layered mate- 
rials. The simplest form for the nonlocal correlation- 
energy part to such a functional is 

E? = - / d 3 rd 3 r' n(r)(j){f,f")n(f'), (1) 



where <p(r, r') is some given, general function depending 
on r — r ' and the densities n in the vicinity of r and r 1 . 
It is approximately derived here and applied to some key 
cases, with results that give promise for broader applica- 
tions. 

We start with the same approximation scheme used for 
layered systems 0, , and divide the correlation energy 
into two pieces, 

E c [n]=E° c [n} + E?[n], (2) 

which are treated in different approximations. In par- 
ticular, as discussed in Ref. 0, we treat the second term 



in the full potential approximation (FPA), which is ex- 
act at long distances between separated fragments, and 
therefore adopt Eq. (25) of Ref. 
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where \ is the density response to a fully self-consistent 
potential with long-range, interfragment spectator || 
contributions omitted |^. V is the interelectronic 
Coulomb interaction, e an appropriately approximated 
dielectric function, and u the imaginary frequency. While 
Eq. © is taken to be the definition of E£ , we will show 
that Eq. can be obtained with suitable approxima- 
tions. The first term E®, defined by Eqs. and ©, is 
also nonlocal; however, with the long-range vdW terms 
treated separately, it seems a reasonable approximation 
to treat E® in the LDA, and this is what we do. There 
is no double counting, because for a uniform system 
1 — Vx = e, which implies that the LDA for i?" 1 van- 
ishes. This is also the key to a seamless theory. 

For layered systems the scheme was made tractable by 
use of the lateral average of the densities to calculate 
the interplanar contribution from Eq. (|3J). For general 
geometries, we make the scheme tractable by expanding 
Eq. J3J to second order in S = 1 — e _1 , obtaining 
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This vanishes in the uniform limit as it must. 

To assess the nature of this approximation, we con- 
sider its prediction for the long-range interaction be- 
tween two fragments. As in Ref. 6, we obtain the 
standard formula for the coefficient of the R~ 6 term 
in terms of the frequency dependent polarizability ten- 
sor o.Pj(uj) for each fragment p = 1,2. We find that 
c?.. = 5 i:j J d 3 x J d 3 y S p (x, y)/47r. We note that S(x, y)/Air 
is the polarization P at point x induced by the D-field 
at point y, i.e., D{y) (in this case the field from fluctua- 
tions on the other fragment) . The integration over x then 
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gives the dipole moment, and the integration over y im- 
plies that this Z?-field is assumed to be spatially constant. 
This means that, although a local-field correction is in- 
cluded, it is not calculated self-consistently. It is the same 
approximation for the electrodynamics that was used to 
calculate the Cg coefficient successfully over a large set 
of atomic and molecular pairs 0, , an approximation 
also independently developed by another group || . Fully 
consistent electrodynamics gave a perceptible, but not 
dramatic, statistical improvement |9j for those atomic 
and molecular dimers tested. However, one should use 
an improved approximation in cases where the fragment 
polarizabilities have significant anisotropy and also for 
the asymptotic interaction between two parallel surfaces, 
which is not given exactly by this approximation. In 
future versions of the functional, it will probably be pos- 
sible to eliminate this ^-expansion, keeping only the ex- 
pansion in terms of the difference between e and 1 — V\- 
The retention of the anisotropic components of e should 
also be possible. 

In order to evaluate Eq. Q), we need a simple approx- 
imation for S, as a functional of the density. This choice 
is constrained by a number of exact relationships. In 
a plane-wave representation, S^?'' one nas the require- 
ments (i) S^iIuj) — > — (4ire 2 /mLL> 2 )n q -- q *i at large fre- 
quencies (the /-sum rule) , where is the Fourier trans- 
form of the density; (ii) du Sq^{iu) — > 8-K 2 Ne 2 /q 2 
for large q, where N is the electron number, to reproduce 
the exactly known self-correlation; (iii) S^^/ = S-^>^^ 
for time reversal invariance; (iv) a finite S^^i{uj) for 
vanishing q or q' at all non-zero values of w, to give 
an exchange-correlation hole with the correct volume 
(charge conservation) . 

An approximate S inspired by the plasmon-pole model 
successfully applied earlier takes S^* = \[S^i + 
d w here 
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We will take 0J q (r) to be a function of the local density 
at point r and its gradient. The above S satisfies all the 
constraints provided that to q — > q 2 /2m for large q. 

To facilitate the numerical evaluation we choose an oj q 
which depends on a single length scale /, and switches 
from its small-q form to q — \/2ml 2 to its large-g form 
(above) when q ~ l/l. We make an arbitrary choice for 
the switching function letting [1 — }ui q = q 2 /2m. 

The quantity I will be a function of position, and for 
convenience later we define q 2 = j/l 2 , where 7 = Att/9. 
Hence, letting h(y) = 1 — e~ ly , we may write 
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It would of course be preferable to use a form with a 
second length scale, as would be provided by a linear 



term oc qvp in the u q dispersion 0, 0, ^} , but this source 
of error is mitigated by a continuously variable choice of 
go based on the density. We make this choice so the 
exchange-correlation energy density e° c (r) defined by 



E° c = d 3 rel(r)n(r) 



(7) 



produced by this choice of qo(r) corresponds to that of a 
full calculation. 

The expression for E® c corresponding to the approxi- 
mation |PJ for E£ l is simply 
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where E se n subtracts off the internal Coulomb self-energy 
of each electron. As was done previously, we expand to 
lowest order in S in the second equality. Substituting for 
S using JSJ), integrating over u = —iu>, and using Eq. (JTJ 
gives 
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where the second term in the brackets is the self-energy 
subtraction written explicitly. Upon substitution from 
©, one finds 

2 f -f\ foe o 2 

4(0 = —/ dy[h(y)-l] = -—q (f). (10) 

Approximations for e xc are conveniently expressed as 
their ratio to the LDA exchange value s^° A — 
—2>e 2 k-p/^ir 1 where kp = in 2 n. Equation (|1(J|) then im- 
plies that the local value of the parameter go is simply 
given by the local value of hp modulated by an easily 
understood energy ratio, that is 
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Equation ijTT|) is used to determine the qo value to be 
used in Eq. ©, continuously as a function of position. 
For this purpose we use LDA with gradient corrections 
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where Z a b = —0.8491. This is the contribution labeled 
"screened exchange" in the table in the review ofRef.El 
where the relationship to the original work of various au- 
thors, who obtained this quantity from first principles, is 
discussed. The remaining contribution Z c , labeled "fluc- 
tuation" in Ref.0, is part of E™ 1 , and hence inappropri- 
ate to include as part of E® c - It comes from a term in 
perturbation theory that gives a long range vdW-like in- 
teraction if taken at long range instead of in an expansion 
in gradients |l2T ]. 
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Writing Eq. Q in a plane-wave representation gives 

This may be straightforwardly, but tediously, expressed 
in the form of Eq. Q , where the kernel <f> is given by 



<p(r, r ') 

where 
T(w,x,y,z) 
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(15) 
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(w + y) (x + z) (w + z)(y + x) 



W(a, b) — 2[(3 — a 2 )&cos&sina + (3 — 6 2 )acosasin6 
+ (a 2 + b 2 — 3) sinasinfr — 3a6cosacos6]/a 3 6 3 . (16) 

The quantities v and i/ are given by v(y) — y 2 /2h(y /d) 
and u'{y) — y 2 /2h(y/d'), with d = \f — r'\q (r) and 
d' = \f— r'\qo(r'), where qo is given by Eq. ifTTl) . The 
kernel thus depends on r and r ' only through d and 
d! , so that can be tabulated in advance in terms of 
these two variables, or better yet in terms of the sum 
and difference variables D and 5 defined by d = D(l + S) 
and d' = D(l - S). Then < D < oo and < |<5| < 1. 
For large d and d' , the asymptotic form is 
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where C — 12(47r/9) 3 me 4 . In Fig. ^ we show a plot of 
4ttD 2 4> vs. -D for several values of S. The integral of the 
5 = curve vanishes as it must. 

The numerical work uses Eq. (J2J for the correlation 
functional coupled with the Zhang- Yang revPBE 01 ex- 
change functional. This choice is motivated by the work 
of Wu et al. ' which pointed out that a more standard 
GGA predicts substantial binding in rare gas dimers from 
exchange alone, a feature absent for exact Hartree-Fock 
exchange. We found that revPBE exchange does not 
have this property, so by using it we insure that vdW 
binding, a correlation effect, actually comes from the cor- 
relation term in our approximation scheme |l5l |. 

In Table [I] our calculated values for the coefficient Cq 
in the asymptotic interaction —Cq/R 6 between the ele- 
ments of several dimers are compared with those from 
previous related calculations as well as reference values, 
which also have some uncertainty. The new values track 
those from Ref. which used the same approximation for 
the electrodynamics, and are only slightly inferior to the 



values obtained in Ref. 0, where self-consistent electro- 
dynamics was used. 

Figure |3 shows the calculated binding-energy curves 
as functions of separation for Ar and Kr dimers. The 
comparison with experimental values for the binding en- 
ergy and distance [lg illustrate the promise of vdW-DF 
for such systems. The binding-energy curves for a ben- 
zene dimcr in the atop-parallel configuration in Fig. 
illustrate the importance of choosing the right GGA fla- 
vor, i.e. revPBE [lflj . to avoid erroneous attraction in 
exchange-only accounts, as discussed in Refs. and 
Figure 0] illustrates the relative agreement between mod- 
ern wavefunction based calculations 0] and our vdW- 
DF method, and the importance of not using the GGA 
alone. 

The above moderate successes of the vdW-DF pro- 
posed here suggest that its use, along with future im- 
provements, may be a way to proceed for calculating 
properties of vdW bound molecules that are too large 
for wavefunction based methods to be useful. 
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FIG. 1: The kernel <f> in Eq. Q as a function of the di- 
mensionless D parameter for several values of the asymmetry 
parameter 8, as defined in the text. 
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FIG. 2: Calculated interaction energy between two Ar atoms 
(dashed curves) and between two Kr atoms (solid curves). 
The experimental equilibrium data [lr3 | are shown for com- 
parison with the full vdW-DF, and the GGA predictions in 
the revPBE flavor are also shown. 
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FIG. 3: Interaction energy between two benzene molecules in 
the atop-parallel configuration as predicted by the exchange- 
only part of two GGA functionals. The full Hartree-Fock 
results are shown for comparison. 
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TABLE I: Ce values for dimers (Rydberg atomic units). 
Present: from Eq. 1171 . Similar: calculations from Ref. .6 
using the same electrodynamics approximation. Unified: cal- 
culations from Ref. using self-consistent electrodynamics. 
Reference: sources cited in Ref. 0- 
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Ne 
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Ar 


124 


126 


143 


134 


Kr 


238 


245 


291 


266 


Xe 


516 


520 


663 


597 


Mg 


1598 


1513 


1230 


1240 
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